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I.  Introduction 

Two  papers  by  Broadhead  et  al.  (1996, 1997)  have  shown  that  the  Wiggins  blind 
deconvolution  algorithm  (Wiggins,  1978;  Walden,  1985)  can  successfully  be  used  to 
estimate  a  source  signature  in  a  realistic  underwater  acoustics  application.  More  recently, 
Broadhead  and  Pflug  (1998)  have  shown  that  the  deconvolution  algorithm  developed  by 
Cabrelli  (1984)  shows  comparable  results  to  the  well-known  Wiggins  algorithm  in  the 
absence  of  noise,  and  superior  results  in  the  presence  of  noise.  Both  of  the  algorithms 
address  the  blind  deconvolution  problem,  i.e.,  they  attempt  to  recover  either  a  source 
signal,  s( t),  or  an  impulse  response  function,  h(t),  from  a  received  signal,  r(t).  The 
received  signal  is  related  to  the  source  signal  and  impulse  response  function  by  the  process 
of  convolution 

r( t)  =  J  s(t)h( r  -  t)dt  , 

which  is  denoted  more  conveniently  as  r  =  s*h.  The  problem  is  said  to  be  "blind"  when 
neither  the  source  nor  the  impulse  response  function  is  known. 

The  Cabrelli  and  Wiggins  algorithms  are  similar  in  that  they  seek  a  sparse  representation 
for  the  impulse  response  function  that  can  be  used  to  achieve  a  source  estimate,  with  the 
sparseness  measured  by  different  norms  for  the  two  methods.  A  sparse  impulse  response 
function  is  characterized  by  a  few  large  amplitudes  interspersed  with  a  large  number  of 
small.  Obviously,  the  success  of  these  methods  depends  on  how  well  the  propagation  path 
for  a  received  signal  satisfies  the  sparseness  criteria.  In  many  applications,  ambient  noise 
is  a  significant  corrupting  factor,  and  it  is  therefore  important  to  consider  how  such  an 
environment  might  affect  a  deconvolution  algorithm,  as  done  in  Broadhead  and  Pflug 
(1998). 

Alternate  norms  in  the  Wiggins  algorithm  have  been  investigated  briefly  by  Broadhead 
et  al.  (1997)  for  underwater  acoustics  applications  and  by  Nandi  et  al.  (1997)  for 
nondestructive  laser  testing  of  materials.  Here,  alternate  norms  are  investigated  for  the 
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more  noise-resistant  Cabrelli  algorithm,  with  results  using  various  norms  for  the  Wiggins 
algorithm  included  for  completeness. 


II.  Brief  Review  of  Methods 

A  received  signal,  x(t),  can  be  written  as  the  convolution  of  a  source  signal  with  an 
impulse  response  function,  x(7)  =  s(t)*h(t).  The  goal  is  to  find  a  filter  that,  upon 
application  to  the  received  signal,  produces  a  good  estimate  of  the  source  signal  or  impulse 
response  function,  i.e.,  find  fit)  such  that  the  impulse  response  estimate,  h{t),  is  given 
by 


h(t)  =  fit )  *  x(t)  =  f{t)  *  [s(r)  *  hit)]. 

If  h(t )  is  a  good  source  estimate,  i.e.,  h(t)*h~\t)  ~  Sit),  then 


h{t)  =  fit)*  [s(t)*h(t)] 

hit )  *  h~x  it)  =  f(t)  *  sit)  *  [hit)  *  h~'  (r)] 

Sit)  =  fit)*  sit), 

and  the  source  estimate  is  sit)  =  f~\t) . 

A  very  brief  review  of  the  Cabrelli  and  Wiggins  blind  deconvolution  methods  are 
included  here.  Although  the  Wiggins  method  requires  an  iterative  solution,  and  the 
Cabrelli  method  does  not,  computational  intensities  required  to  achieve  solutions  are 
similar. 


A.  Cabrelli  Algorithm 

The  Cabrelli  algorithm  uses  a  measure  of  sparseness  or  simplicity  called  the  D-norm, 
defined  by 


Diy)  =  max 


for  a  real-valued  vector  y  of  length  m,  where 
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M=  lyl 


op 


V  t 


is  the  Euclidean  norm,  which  is  related  to  the  variance.  Here,  y .  =  f.  *  x  is  a  vector  of 
length  m  =  n  + 1  —  1,  xis  the  input  signal  of  length  n,  and  f(  is  a  filter  of  fixed  length  /  (to 


be  determined).  The  blind  deconvolution  problem  is  solved  by  maximizing  D(y)  over  all 
nonzero  filters  f(  ,  i.e.,  maximizing  the  sparseness  of  the  impulse  response  function. 
Mathematically,  D( y)  is  maximized  by  differentiating  D{ y)  with  respect  to  f ,  for  each  i 
and  equating  to  zero  (finding  the  extremal  point  of  D{ y)).  Differentiation  leads  to  the 
matrix  formulation 


f 


1r  •  f  =  x' 


where  x1  -  xj+v  ...,  xj+lA  j  is  an /-length  subset  of  the  n -length  input  channel  x,  R 
is  the  Toeplitz  autocorrelation  matrix  defined  by 


R  = 


r0  ri 

r,  r0 


VrM 


LM 


L 1-2 


ri  f0  J 


and  f .  is  the  set  of  filters  found  by  calculating  f,  =  R1xi .  Once  the  filters  are  found,  they 
are  convolved  with  the  input  ( y(  =  f  ,*x )  to  generate  m  =  n  + 1- 1  potential  impulse 
response  function  estimates.  The  i  for  which  the  D-norm  of  y .  is  maximum  indicates  that 
the  impulse  response  estimate  and  the  corresponding  source  estimate  are  h(t)  =  y,  (r),  and 
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B.  Wiggins  Algorithm 

As  in  the  Cabrelli  algorithm,  the  goal  of  the  Wiggins  algorithm  is  to  find  the  filter  f ., 

that,  upon  convolution  with  the  input  signal  x,  produces  an  estimate  of  the  impulse 
response  function  and  source  signal  (h(t)  =  f(t)*x(t)  and  s(t)  =  f~\t) ).  Instead  of  the 
D-norm,  Wiggins  uses  the  Varimax  norm,  or  V-norm,  as  a  measure  of  sparseness, 


Differentiating  V(y)  with  respect  to  the  filter  coefficients  and  equating  to  zero,  a  set  of 

equations  is  obtained  which  can  be  rewritten  in  matrix  form  as 

R(f)-f  =  g(f) 


where  R(f )  is  the  Toeplitz  autocorrelation  matrix,  and  g(f )  is  a  column  vector 
proportional  to  y3*x.  Choosing  an  initial  filter  f°  =(0,..., 0,1,0, ...,0),  an  iterative 
algorithm  can  be  generated  by  taking 

r+,={R(r)}-,g(r). 

A  convergence  criteria  on  the  estimated  source  or  filter  is  used  to  terminate  the  iterations. 


III.  n-th  Order  Norms  for  Deconvolution 

It  is  shown  in  Broadhead  and  Pflug  (1998)  and  Nandi  et  al.  (1997)  that  the  Wiggins 
algorithm  can  easily  be  extended  to  used  alternate  norms,  i.e.,  the  n-th  order  normalized 
moment  instead  of  the  fourth-order  normalized  moment,  or  V-norm.  The  Cabrelli 
algorithm  can  similarly  be  modified  by  simply  substituting  alternate  norms  for  the  D-norm 
in  the  implementation  of  the  originally-derived  formulation.  Unlike  the  formally-derived 
extension  of  the  Wiggins  algorithm  included  in  Broadhead  and  Pflug,  the  extension  of  the 
Cabrelli  algorithm  to  n-th  order  norms  is  heuristic,  based  on  the  observation  that  the  D- 
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norm  tends  to  echo  the  signal  kurtosis,  and  that  therefore  norms  based  on  moments  (fourth 
and  otherwise)  might  be  useful. 

IV.  Simulated  Data 

Two  source  signals  and  two  impulse  response  functions  are  used  to  compare  the 
performance  of  n-th  order  norms  in  the  Cabrelli  and  Wiggins  deconvolution  algorithms. 
The  first  source  (sj)  is  a  pulse-type  signal  and  the  second  source  (s'2)  is  an  exponentially- 
damped  sinusoid,  as  shown  in  Fig.  1.  The  first  impulse  response  function  (hf)  is  a  series 
of  five  positive  spikes  with  overall  skewness  equal  to  1 1.4  and  kurtosis  equal  to  136.5. 
Since  a  skewed  signal  is  by  definition  also  nonGaussian,  the  first  impulse  response  has  by 
necessity  both  significant  skew  and  kurtosis.  The  second  impulse  response  function  (h.2) 
is  a  series  of  five  positive  and  negative  spikes  with  overall  skewness  equal  to  0.9  and 
kurtosis  equal  to  81.8.  That  is,  /z2  is  fairly  symmetric,  but  nonGaussian.  The  second 
impulse  response  is  representative  of  multipath  propagation  that  might  occur  in  an 
underwater  acoustic  waveguide  with  air  at  the  surface,  for  example.  In  contrast,  the  first 
might  represent  propagation  in  an  ocean  environment  with  a  smooth  ice  cover.  These  two 
impulse  response  functions  provide  an  opportunity  to  judge  the  relationship  between  the 
normalized  moments  of  an  impulse  response  function  and  the  order  of  the  norm  used  in  the 
deconvolution  algorithms. 

Each  signal  is  convolved  with  each  impulse  response  function  to  create  a  set  of 
simulated  received  signals  for  input  into  the  deconvolution  algorithms.  These  input  signals 
are  shown  in  Fig.  2.  The  spikes  in  hj  and  h.2  are  spaced  closely  enough  that  multipath 
arrivals  of  the  signal,  especially  S2,  are  not  spatially  resolved. 
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Figure  2.  Input  signals  for  the  deconvolution  algorithms,  where  sj  and  si  denote  the 
first  and  second  source  signals,  hi  and  hi  denote  the  first  and  second  impulse  response 
functions,  and  denotes  the  convolution  operation.  The  four  input  signals  are  (a)  sl*hl, 
(b)  5,  *  /zj ,  (c)  s2*hi,  and  (d)  s2*h7. 


V.  Algorithm  Parameters  and  Performance  Evaluation 

For  the  Cabrelli  and  Wiggins  algorithms,  the  3rd,  4th,  5th,  and  6th  order  normalized 
moments  are  evaluated,  with  the  normalized  moments  defined  by 

2>-50" 

m  =  - 

*  mn 

for  a  signal  sampled  with  N  points  and  mean  y.  For  the  Wiggins  algorithm  with  n  =  4, 
this  is  essentially  equivalent  to  the  V-norm,  and  is  referred  to  as  such. 

Filter  lengths  from  1  to  50  are  tested.  For  most  of  the  cases  tested  here,  at  least  one  of 
these  filter  lengths  is  sufficient  to  produce  a  good  source  estimate.  In  practice,  the  filter 
length  is  unknown,  but  a  set  of  filter  lengths  could  routinely  be  tested  to  produce  a  set  of 
possible  solutions,  which  would  then  have  to  be  evaluated  in  some  systematic  manner. 

To  evaluate  algorithm  performance,  a  rather  simple,  but  effective,  measure  is  used.  The 
correlation  coefficient  (cc)  between  the  source  estimate,  s(t )  and  the  true  source,  s(t) 
given  by 

max 

(IXMArfO-’MArf  ' 

This  quantity  is  bounded  between  zero  and  one,  with  a  value  of  one  indicating  that  the 
source  estimate  is  equal  to  the  true  source. 

The  Wiggins  algorithm  requires  a  convergence  criteria  for  the  iterative  solution  to  the 
nonlinear  system  of  equations.  This  is  chosen  to  be  either  the  point  at  which  the 
correlation  coefficient  between  the  current  and  previous  source  estimate  is  0.9999,  or  at 
100  iterations.  When  the  goal  is  to  estimate  the  impulse  response  function  rather  than  the 
source  signal,  the  criteria  may  more  appropriately  be  placed  on  the  estimated  impulse 
response  function.  Although  prewhitening  is  sometimes  required  for  inversion  of  the 
autocorrelation  matrix,  it  was  not  needed  in  these  simulations.  However,  restricting  the 
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input  and  output  signals  to  an  estimate  of  the  source  signal  passband  was  required  in  one 
case  to  achieve  successful  results  with  the  Wiggins  algorithm.  This  information  would 
generally  be  available  for  high  signal-to-noise  ratio  signals,  and  is  not  overly  restrictive. 

VI.  Deconvolution  Results 

This  section  contains  the  results  of  the  Cabrelli  and  Wiggins  algorithm  performance. 
No  noise  is  included  in  the  simulations,  although,  as  shown  in  Broadhead  and  Pflug 
(1998)  and  claimed  in  Cabrelli  (1984),  the  Cabrelli  method  appears  to  be  more  robust  to 
additive  noise  than  the  Wiggins  algorithm. 

The  results  are  shown  in  two  forms.  The  first  form  is  a  figure  depicting  the  correlation 
coefficient  between  the  source  estimate  and  the  true  source  at  each  filter  length.  It  is 
generally  desirable  that  good  source  estimates  be  produced  for  many  filter  lengths,  as  the 
required  filter  length  will  be  unknown  in  practice.  The  second  form  of  results  is  a  figure 
depicting  the  best  source  estimate  over  filter  length  for  each  method  and  norm.  This  allows 
a  visual  evaluation  of  the  results  and  provides  a  guideline  for  determining  the  significance 
of  the  correlation  coefficients.  For  comparison,  the  correlation  coefficients  between  the 
unprocessed  input  signals  and  the  true  sources  are  given  in  the  Table. 


Input  Signal 

Correlation  Coefficient 

0.8834 

0.7247 

S2*\ 

0.9466 

S2**h 

0.7961 

Table.  Correlation  Coefficients  for  the  Original  Input  Signals  (Before  Deconvolution) 

A.  Signal  1 

The  correlation  coefficients  versus  filter  length  for  the  Cabrelli  algorithm  are  given  in 
Fig.  3.  In  this  and  similar  figures,  the  correlation  coefficients  between  the  unprocessed 
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signals  and  the  true  sources  from  the  Table  are  depicted  by  horizontal  dashed  black  lines. 
When  the  first  impulse  response  is  used  to  create  a  received  signal  for  input  into  the 
deconvolution  routine  (s,*/z,),  the  D-norm,  3rd,  4th,  5th,  and  6th  order  norms  produce 
good  source  estimates  for  many  filter  lengths.  In  this  case,  good  solutions  are  those  with 
correlation  coefficients  of  0.8960  and  higher,  and  are  only  slightly  better  than  the 
correlation  coefficient  for  the  unprocessed  signal,  which  is  0.8834.  Only  the  3rd  moment 
norm  produces  a  significantly  improved  source  estimate  with  cc  =  0.9496  (see  Fig.  4),  and 
only  for  one  filter  length.  The  Wiggins  algorithm  with  each  norm  produces  at  least  one 
source  estimate  that  is  better  than  the  input  signal,  as  shown  by  the  results  in  Figs.  5  and  6, 
with  the  best  estimate  given  by  the  5th  order  norm,  cc  =  0.9506.  Most  of  the  good 
solutions  are  only  slight  improvements,  however,  and  they  occur  at  fewer  filter  lengths  in 
the  Wiggins  results  compared  to  the  Cabrelli  results. 

For  the  second  impulse  response,  ( sx  *  h^),  the  Cabrelli  and  Wiggins  algorithms 
produce  improved  source  estimates  over  several  filter  lengths.  In  both  algorithms,  the  5th 
order  norm  gives  the  best  solution  with  cc  =  0.8941  for  the  Cabrelli  algorithm  and 
cc  =  0.8827  for  the  Wiggins  algorithm.  These  are  significant  improvements  over  the 
unprocessed  input  signal  correlation  coefficient  of  0.7247.  Note,  however,  that  the  D- 
norm  and  the  Cabrelli  algorithm  gives  produces  improved  source  estimates  at  more  filter 
lengths  than  the  V-norm  and  the  Wiggins  algorithm,  and  the  best  source  estimate  from  the 
D-norm  has  cc  =  0.8706,  which  is  comparable  to  the  best  source  estimate  from  the  V-norm 
with  cc  =  0.8787. 

B.  Signal  2 

Results  for  the  second  test  signal  are  shown  in  Figs.  7  and  8  for  the  Cabrelli  algorithm 
and  in  Figs.  9  and  10  for  the  Wiggins  algorithm.  The  s2*hl  input  signal  is  very  similar  to 
the  true  source,  having  cc  =  0.9466.  Even  so,  both  the  Cabrelli  and  Wiggins  algorithms 
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with  all  order  norms  produce  many  source  estimates  over  the  filter  lengths.  All  order 
norms  in  the  Cabrelli  algorithm  produce  best  source  estimates  with  cc  =  0.9971,  and  all  the 
best  source  estimates  are  above  cc  =  0.9963  for  the  Wiggins  algorithm. 

The  more  significant  distortions  in  the  inputs  signal  does  not  prevent  either 
algorithm  from  producing  excellent  source  estimates,  but  they  occur  over  only  about  one- 
half  of  the  filter  lengths  tested.  Also,  for  the  Cabrelli  algorithm,  only  the  D-norm,  4th  and 
6th  moment  norms  produces  significantly  improved  source  estimates,  with  cc  =  0.9951  for 
each,  compared  to  the  input  signal  correlation  coefficient  of  0.7961.  The  3rd  and  5th 
moment  norms  produce  source  estimates  that  appear  to  be  reversed.  Since  there  is  no 
obvious  reason  why  this  should  occur  in  either  the  Cabrelli  or  Wiggins  algorithms,  these 
occurrences  are  assumed  at  this  point  to  be  arbitrary.  However,  we  note  that  if  these 
source  estimates  were  reversed,  the  correlation  coefficients  would  increase,  from  0.7855  to 
0.9322  for  the  3rd  order  norm  and  from  0.7814  to  0.9586  for  the  5th  order  norm. 

Whether  this  is  likely  to  occur  for  other  signal  types  is  as  yet  unknown.  In  contrast  to  the 
Cabrelli  results,  all  the  norms  in  the  Wiggins  algorithm  produce  good  source  estimates, 
and  have  correlation  coefficients  greater  than  0.9966. 
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Filter  Length 


Figure  3.  Correlation  coefficients  between  the  estimated  pulse  source,  sx,  and  the  true 
source  versus  filter  length  using  the  Cabrelli  algorithm  for  input  signals  (a)  s^l\,  and  (b) 

s*hi. 


12 


Source  1  -  Cabrelli 


D-norm 


D-norm 


Sample  Number 


Sample  Number 


Figure  4.  Best  estimated  source  signals  for  input  signals  (a)-(e)  ,  and  (f)-(j)  s^h^, 

resulting  from  the  Cabrelli  algorithm  with  various  norms.  The  correlation  coefficient  (cc) 
between  the  source  estimate  and  the  true  source  is  given  within  the  plot. 
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Filter  Length 


Figure  5.  Correlation  coefficients  between  the  estimated  pulse  source,  sx  ,  and  the  true 
source  versus  filter  length  using  the  Wiggins  algorithm  for  input  signals  (a)  s^*l\ ,  and  (b) 
s^h,. 
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Figure  6.  Best  estimated  source  signals  for  input  signals  (a)-(d)  ,  and  (e)-(h) 

resulting  from  the  Wiggins  algorithm  with  various  norms.  The  correlation 
coefficient  (cc)  between  the  source  estimate  and  the  true  source  is  given  within  the  plot. 
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Figure  7.  Correlation  coefficients  between  the  estimated  damped  sinusoid  source,  s2 
and  the  true  source  versus  filter  length  using  the  Cabrelli  algorithm  for  input  signals  (a) 
and  (b) 
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Figure  8.  Best  estimated  source  signals  for  input  signals  (a)-(e)  s2*h^,  and  (f)-(j)  s2*h2, 
resulting  from  the  Cabrelli  algorithm  with  various  norms.  The  correlation  coefficient  (cc) 
between  the  source  estimate  and  the  true  source  is  given  within  the  plot. 
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Figure  9.  Correlation  coefficients  between  the  estimated  damped  sinusoid  source,  s2 , 
and  the  true  source  versus  filter  length  using  the  Wiggins  algorithm  for  input  signals  (a) 
s2*I\,  and  (b)  s^h,. 
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Figure  10.  Best  estimated  source  signals  for  input  signals  (a)-(d)  ,  and  (e)-(h) 

si*K'  resulting  from  the  Wiggins  algorithm  with  various  norms.  The  correlation 
coefficient  (cc)  between  the  source  estimate  and  the  true  source  is  given  within  the  plot. 
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VI.  Conclusions 


From  these  results,  it  appears  that  there  is  generally  little  justification  for  using  norms 
other  than  the  D-norm  for  the  Cabrelli  algorithm,  or  the  V-norm  for  the  Wiggins  algorithm. 
While  the  results  using  other  norms  are  better  in  some  cases,  meaning  that  the  best  source 
estimate  is  more  similar  to  the  true  source,  or  good  source  estimates  are  produced  more 
consistently  with  varying  filter  length,  no  predictable  pattern  emerges  to  provide  guidelines 
as  to  when  which  norm  will  work  best  in  most  cases.  One  could  routinely  use  more  than 
one  norm  for  estimation,  but  the  number  of  possible  solutions  then  increases  significantly, 
and  unless  one  has  an  effective  method  of  handling  the  large  number  of  solutions  generated 
by  multiple  norms  and  multiple  filter  lengths,  there  may  be  no  improvement.  The 
exception  is  an  application  in  which  the  impulse  response  function  is  known  to  consist  of 
all  positive  spikes,  and  the  source  is  similar  to  the  smooth,  symmetric  wavelet-type  source 
( 5, )  (what  constitutes  "similar"  has  not  been  determined).  In  this  case,  the  3rd  and  5th 
order  norms  appear  to  work  better  in  the  Cabrelli  algorithm  than  the  D-norm  or  even  order 
norms.  Note  that,  excepting  the  cases  in  which  the  source  estimate  appears  time-reversed, 
each  norm  and  each  algorithm  produced  a  source  estimate  superior  to  the  original  input 
signal. 
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